arXiv: 1507.08612v2 [stat.ME] 31Jul2015 


Likelihood-free inference in high-dimensional models 


Athanasios Kousathanas 1,2 , Christoph Leuenberger 3 , Jonas Heifer 4 , Mathieu 
Quinodoz 5 , Matthieu Foil 5 ' 6 , and Daniel Wegmann 1,2 

1 Department of Biology and Biochemistry, University of Fribourg, Fribourg, 

Switzerland 

2 Swiss Institute of Bio informatics, Lausanne, Switzerland 
3 Department of Mathematics, University of Fribourg, Fribourg, Switzerland 
4 Massachusetts Institute of Technology (MIT),Cambridge MA, USA 
5 School of Life Sciences, Ecole Polytechnique Federate de Lausanne (EPFL), 

Lausanne, Switzerland 

international Agency for Research on Cancer, Lyon, France 


Abstract 

Methods that bypass analytical evaluations of the likelihood function have be¬ 
come an indispensable tool for statistical inference in many fields of science. These 
so-called likelihood-free methods rely on accepting and rejecting simulations based 
on summary statistics, which limits them to low dimensional models for which the 
absolute likelihood is large enough to result in manageable acceptance rates. To get 
around these issues, we introduce a novel, likelihood-free Markov-Chain Monte Carlo 
(MCMC) method combining two key innovations: updating only one parameter per 
iteration and accepting or rejecting this update based on subsets of statistics suffi¬ 
cient for this parameter. This increases acceptance rates dramatically, rendering this 
approach suitable even for models of very high dimensionality. We further derive that 
for linear models, a one dimensional combination of statistics per parameter is suffi¬ 
cient and can be found empirically with simulations. Finally, we demonstrate that our 
method readily scales to models of very high dimensionality using both toy models 
as well as by jointly inferring the effective population size, the distribution of fitness 
effects of new mutations (DFE) and selection coefficients for each locus from data of 
a recent experiment on the evolution of drug-resistance in Influenza. 

Significance Statement 

The goal of statistical inference is to learn about the parameters of a model that 
led to the data observed. In complex models, this is often difficult due to a lack of 
analytical solutions. A popular solution is to replace direct calculations with computer 
simulations, but the inefficiency of sampling algorithms currently restricts this to low¬ 
dimensional models. Here we construct a novel approach that exploits the observation 
that the information about a parameter is often contained in a subset of the data. This 
approach readily scales to high dimensions and enables inference in more complex 
and possibly more realistic models. It allowed us, for instance, to accurately infer 
parameters underlying the evolution of drug-resistance in Influenza viruses. 
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The past decade has seen a rise in the application of Bayesian inference algorithms that 
bypass likelihood calculations with simulations. Indeed, these generally termed likelihood- 
free or Approximate Bayesian Computation (ABC) [T| methods have a wide range of ap¬ 
plications ranging from communications engineering to population genetics [2]. This is 
because many scientific fields employ complex models for which likelihood calculations are 
intractable, thus necessitating inference through simulations. 

Let us consider a model Ai that depends on n parameters 9 , creates data D and has 
the posterior distribution 


ir(9\D) 


F(D\9)tt(9) 

I P(D\9)tt(9) ’ 


where tt(9) is the prior and P(D|0) is the likelihood function. ABC methods bypass 
the evaluation of P(Z)|0) by performing simulations with parameter values sampled from 
7 t(9) that generate D, which in turn is summarized by a set of m-dimensional statistics s. 
The posterior distribution is then evaluated by accepting such simulations that reproduce 
the statistics calculated from the observed data (s 0 bs) 


7r(<?|s) 


P(s = S obs \9)TT(6) 

/P(s = s obs \9)n(9)' 


However for models with m » 1 the condition s = s obs might be too restrictive and 
requiring a prohibitively large simulation effort. Therefore, an approximation step can be 
employed by relaxing the condition s = s obs to || s — s obs ||< 6 , where || x — y || is an 
arbitrary distance metric between x and y and 5 is a chosen distance (tolerance) below 
which simulations are accepted. The posterior 7r(0|s) is thus approximated by 


_/0i,x = gCI s ~ s obs ||< S\6)tt(6) 

/P(|| s - s obs ||< d^irte) 1 

An important advance in ABC inference for models of low to moderate dimension was 
the development of methods coupling ABC with Markov Chain Monte Carlo (MCMC) 
[3].These methods allow efficient sampling of the parameter space in regions of high like¬ 
lihood, thus requiring less simulations to obtain posterior estimates [3j. The original 
ABC-MCMC algorithm proposed by Marjoram et al. (2003) [3] is: 


1. If now at 9 propose to move to 9' according to the transition kernel q(9'\9). 

2. Simulate D using model Ai with 9' and calculate summary statistics s for D. 

3. If || s — s 0 fe s || < S go to step 4 otherwise go to step 1. 

4. Calculate the Metropolis-Hastings ratio 


h = h(9, 9') = min 



A0')<1{9\9') \ 
n(9)q(9'\9) ) ’ 


5. Accept 9' with probability h otherwise stay at 9. Go to step 1. 


The sampling success of ABC-MCMC is given by the absolute likelihood values, which are 
often very low even for relative large tolerance values 5 . In such situations, the condition 
|| s — s obs ||< 5 will impose a quite rough approximation to the posterior. As a result, the 
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utility of ABC-MCMC is limited to models of relatively low dimensionality (typically up to 10 
parameters; 13 El)- The same limitation applies to the more recently developed sequential 
Monte Carlo sampling methods HUH]. 

To this end, three approaches have been suggested to address models of higher dimen¬ 
sionality with ABC. The first approach proposes an expectation propagation approximation 
to factorize the data space |9|, which is an efficient solution for situations with high dimen¬ 
sional data, but does not directly address the issue of high dimensional parameter spaces. 
The second approach formulates the problem using hierarchical models and proposes to 
first estimate the hyper parameters, and then fixing them when inferring parameters of 
lower hierarchies individually [10|. The third approach consists of first inferring marginal 
posterior distributions on low dimensional subsets of the parameter space (either one |11| 
or two dimensional |12j). and then reconstructing the joint posterior distribution from 
those. The latter two approaches benefit from the lower dimensionality of the statistics 
space when considering subsets of the parameters individually, and hence render the ac¬ 
ceptance criterion meaningful. However, they will not recover the true joint distribution if 
parameters are correlated, which is a common feature of complex models. 


Efficient ABC in high-dimensional models 


Here we introduce a new ABC algorithm that exploits the reduction of dimensionality of 
the summary statistics when focusing on subsets of parameters, but couples the parameter 
updates in an MCMC framework. As we prove below, this coupling ensures that our 
algorithm converges to the true joint posterior distribution even for models of very high 
dimensions. 

Let us define the random variable T, = T,(s) as an mi-dimensional function of s. We 
call Tj sufficient for the parameter 0 ? ; if the conditional distribution of s given T; does not 
depend on 6j. More precisely, let ti fibs = Tj(s obs ). Then 


P(^ 8 obs \Tj tj obsi 


(® — 8 obs: Tj — ti,obs {O') 
P(T , i = t i:Obs \0) 

*(8 = S obs \6) 


P(Tj = t ifibs \0) 


— • 9^y8 0 bsi 


where = (#i,..., 0j_i, #i+i,..., 6 n ) is 6 with the i-th component omitted. 


If sufficient statistics Tj can be found for each parameter 0, and their dimension m* 
is substantially smaller than the dimension m of s, then the ABC-MCMC algorithm can be 
greatly improved with the following algorithm which we denoted ABC with Parameter 
Specific Statistics or ABC-PaSS onwards: 

The algorithm starts at time t = 1 and at some initial parameter value 0^\ 

1. Choose an index i = 1 ,... ,n according to a probability distribution (pi,... ,p n ) with 
YhPi = 1 and a ll Pi > 0. 

2. At 6 = 0 (i) propose O' according to the transition kernel qj[0'\0 ) where O' differs 
from 0 only in the i-th component: 

O' = (0i, . . . ,0j_i,0',0 i+ i, . . . , 0 n ). 

3. Simulate D using model Ai with O' and calculate summary statistics s for D. Cal¬ 
culate ti — Ti(s) and ti, obs — Tj(s ob ^ j. 
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4. Let Si be the tolerance for parameter d L . If \\t t — tj )0 b s || < Si go to step 5 otherwise 
go to step 1. 

5. Calculate the Metropolis-Hastings ratio 



6. Accept 6' with probability h otherwise stay at G. 

7. Increase t by one, save a new parameter value 6 ^ = 6 and continue at step 1. 


Convergence of the MCMC chain is guaranteed by 

Theorem 1. Fori = l..n, if Si = 0 andT{ is sufficient for parameter 6i then the stationary 
distribution of the Markov chain is 7r(0|s = s 0 b s ). 

The proof for Theorem [l] is provided in the Appendix. 

The same types of improvements that have been proposed for ABC-MCMC can also be 
used for ABC-PaSS. For example in order to increase acceptance rate we can relax the 
assumption Si = 0 to <5* > 0 and the stationary distribution of the Markov chain is then 
7r(0|||£'-ii|| i < Si, i = 1,..., n) and approximates the posterior distribution 7r(0|s = s Q b s ). 
We can also perform an initial calibration ABC step to find an optimal starting position 
tolerance Si and to adjust the proposal kernel for each parameter [4]. 

Toy model 1: Normal distribution 

We first compared the performance of ABC-PaSS and ABC-MCMC under a simple model: the 
normal distribution with parameters mean (/x) and variance (cr 2 ). Given a sample of size 
n, the sample mean (x) is a sufficient statistic for /x, while both x and the sample variance 
(5 2 ) are sufficient for a 2 [13] . For ABC-MCMC, we used both x and S 2 as statistics. For 
ABC-PaSS, we used only x when updating fi and both x and S 2 when updating a 2 . 

We then compared the accuracy between the two algorithms by calculating the total 
variation distance between the inferred and the true posteriors ( L\ distance). We computed 
L\ under a wide range of tolerances in order to find the tolerance for which each algorithm 
had the best performance (i.e., minimum L\). As shown in Figure [I] panels A and C, 
ABC-PaSS produced a more accurate estimation for fj, than ABC-MCMC. The two algorithms 
had similar performance when estimating a 2 (Figure [T| B and D). 

The normal distribution toy model, although simple, is quite illustrative of the nature 
of the improvement in performance by using ABC-PaSS over ABC-MCMC. Indeed, our results 
demonstrate that the slight reduction of the summary statistics space by ignoring a single 
uninformative statistic when updating /x already results in a noticeable improvement in 
estimation accuracy. This improvement would not be possible to attain with classic di¬ 
mension reduction techniques, such as partial least squares (PLS) since the information 
contained in x and S 2 is irreducible under ABC-MCMC. 

Sufficient Statistics in high-dimensional models 

Decreasing the dimensionality of statistics space is crucial for ABC-PaSS, since we would 
expect an improvement over ABC-MCMC only if we can find per-parameter sufficient statistics 
of a lower dimension than the total number of parameters. However, choosing summary 


4 



statistics is not trivial in any ABC application, as too few statistics are insufficient to 
summarize the data while too many statistics can create an excessively large statistics space 
that worsens the approximation of the posterior PHQU- Therefore, several strategies 
have been developed to reduce the dimensionality of the statistics space |15] , For instance, 
Fearnhead and Prangle |6] suggested a method where an initial set of simulations is used 
to fit a linear model that expresses each parameter 6i as a function of s. These functions 
are then used as statistics in subsequent ABC analysis. 

Here we will adopt Fearnhead and Prangle’s approach to reduce the dimensionality of 
statistics space to a single combination of statistics per parameter. This choice is moti¬ 
vated by the finding that for a general linear model (GLM), a single linear function is a 
sufficient statistic for each associated parameter, as we prove in the following. 

Suppose that, given the parameters 6 , the distribution of the statistics vector s is 
multivariate normal according to the general linear model (GLM) 

s = c + CO + e 

where e ~ A/ r (0, E s ) and for any m x n-matrix C. If the prior distribution of the parameter 
vector is 6 ~ N(Oq, Eg) then the posterior distribution of 0 given s 0 b s is 

0\s 0 bs ~ Af{Dd, D) (1) 

with D = (C / S“ 1 C + Eg 1 ) _1 and d = C / 5]J 1 (s 0 b s — c) + Eg 1 0 q (see e.g. [Ej). We have 
the 

Theorem 2. Let Ci be the i-th column of C and (3i = Ej 1 ^. Moreover, let 

Ti = Ti(s ) = /3'S. 

Then Ti is sufficient for the parameter 6i and the collection of statistics 

T = (Ti, . . . ,T n )' 


yields the same posterior 0 as s. 

The proof for Theorem [2] is provided in the Appendix. In practice, the design matrix 
C is unknown. We can perform an initial set of simulations from which we can infer that: 


Cov(s, Of) = Var(0j)cj. 

- / 

A reasonable estimator for the sufficient statistic r* is then Tj = f3 t s with 

Pi = ( 2 ) 

where and E s g. for i = 1,... ,n are the estimated covariances. 


Toy model 2: General Linear Model (GLM) 

We next compared the performance of ABC-MCMC and ABC-PaSS under GLM models of 
increasing dimensionality n. For all models, we constructed the design matrix C such 
that all statistics are informative for all parameters, while retaining the total information 
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on the individual parameters regardless of dimensionality (see methods). For ABC-MCMC, 
we used all statistics s, while for ABC-PaSS, we employed Theorem [2] and used a single 
linear combination of statistics r, per parameter 0j. As above, we assessed performance 
of ABC-MCMC and ABC-PaSS by calculating the total variation distance {L\) between the 
inferred and the true posterior distribution. We calculated L\ for several tolerances in 
order to find the tolerance where L\ was minimal for each algorithm (see Figure for 
examples with n = 2 and n = 4). Since in ABC-MCMC distances are calculated in the multi¬ 
dimensional statistics space, the optimal tolerance increased with higher dimensionality. 
This is not the case for ABC-PaSS, because distances are always calculated in one dimension 
only (Figure |2]A). 

We found that ABC-MCMC performance was good for low n, but worsened rapidly with 
increasing number of parameters, as expected from the corresponding increase in the di¬ 
mensionality of statistics space (Figure [2^3). For a GLM with 32 parameters, approximate 
posteriors obtained with ABC-MCMC differed only little from the prior (Figure [2^3). In con¬ 
trast, performance of ABC-PaSS was unaffected by dimensionality and was better than that 
of ABC-MCMC even in low dimensions (Figure [2^3). These results support that by consid¬ 
ering low dimensional parameter-specific summary statistics under our framework, ABC 
inference remains feasible even under models of very high dimensionality, for which current 
ABC algorithms are not capable of producing meaningful estimates. 


High-dimensional population genetics inference of natural se¬ 
lection and demography 

One of the major research problems in modern population genetics is the inference of 
natural selection and demographic history, ideally jointly [ T71 fl8] . One way to gain insight 
into these processes is by investigating how they affect allele frequency trajectories through 
time in populations, for instance under experimental evolution. Several methods have 
thus been developed to analyze allele trajectory data in order to infer both locus-specific 
selection coefficients (s) and the effective population size (N e ). The modeling framework 
of these methods assumes Wright-Fisher (WF) population dynamics in a hidden Markov 
setting to calculate the likelihood of the observed allele trajectories give parameters N e and 
s [l9l HJj. In this setting, likelihood calculations are feasible, but very time-consuming, 
especially when considering many loci at the genome-wide scale [21j |. 

To speed-up calculations, Foil et al. m developed an ABC method (WF-ABC), adopting 
the hierarchical ABC framework of Bazin et al. [TO] . Specifically, WF-ABC first estimates 
N e based on statistics that are functions of all loci, and then infers s for each locus indi¬ 
vidually under the inferred value of N e . While WF-ABC easily scales to genome-wide data, 
it suffers from the unrealistic assumption of complete neutrality when inferring N e , which 
is potentially leading to biases in the inference. 

Here we employed ABC-PaSS to infer both N e and locus-specific selection coefficients 
jointly. To reduce the dimensionality of the statistics space, we fitted parameter-specific 
linear combinations as described above, motivated by the frequent and successful applica¬ 
tion of linear approximations in ABC [22j 4]- In addition, we first applied a multivariate 
Box-Cox transformation [23] to increase linearity between statistics and parameters, as sug¬ 
gested by [3], and then assessed the assumption of linearity empirically (Supplementary 
Figure SI). 
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Perfomance of ABC-PaSS in inferring selection and demography 

To examine the performance of ABC-PaSS under the WF model, we inferred N e and s on 
sets of 100 loci simulated with varying selection coefficients. We evaluated the accuracy of 
estimation by comparing the estimated versus the true values of the parameters over 25 
replicate simulations. As shown in Figure [3]A, N e was estimated well over the whole range 
of values tested. Estimates for s were on average unbiased and accuracy was, as expected, 
higher for larger N e (Figure [3^3). Note that since the prior on s was C/[0,1], these results 
imply that our approach estimates N e with high accuracy even when the majority of the 
simulated loci are under strong selection (90% of loci had N e s > 10). Hence, our method 
allows to relax the assumption of neutrality on most of the loci, which was necessary in 
previous studies (|2Tj). 

We next introduced hyper parameters for the distribution of selection coefficients (the 
so called distribution of fitness effects or DFE). Such hyper parameters are computationally 
cheap to estimate under our framework, as their updates can be done analytically and do 
not require simulations. Following |M1122], we assumed that the distribution of the locus- 
specific s is realistically described by a truncated Generalized Pareto Distribution (GPD) 
with location /j = 0 and parameters shape a and scale y (Supplementary Figure S2). 

We first evaluated the accuracy of estimating y and a when fixing the value of the other 
parameter and found that both parameters are well estimated under these conditions (Fig¬ 
ure [3} C and D, respectively). Since the truncated GPD of multiple combinations of y 
and a is very similar, these parameters are not always identifiable. This renders the accu¬ 
rate joint estimation of both parameters difficult (Supplementary Figure 3B-C). However, 
despite the reduced accuracy on the individual parameters, we found the overall shape of 
the GPD to be well recovered (Supplementary Figure 3D-F). Also, N e was estimated with 
high accuracy for all combinations of y and a (Supplementary Figure 3A). 

Analysis of Infuenza data 

We applied our approach to data from a previous study [26j where cultured canine kidney 
cells infected with the Influenza virus were subjected to serial transfers for several genera¬ 
tions. In one experiment, the cells were treated with the drug Oseltamivir, and in a control 
experiment they were not treated with the drug. To obtain allele frequency trajectories 
of all sites of the Infuenza virus genome (13.5 Kbp), samples were taken and sequenced 
every 13 generations with pooled population sequencing. The aim of our application was 
to identify which viral mutations rose in frequency during the experiment due to natural 
selection and which due to drift and to investigate the shape of the DFE for the control 
and drug-treated viral populations. 

Following [26] , we filtered the raw data to contain loci for which sufficient data was 
available to calculate the summary statiatics considered here (see Methods). There were 86 
and 42 such loci for the control and drug-treated experiment, respectively (Supplementary 
Figure S4). 

We then employed ABC-PaSS to estimate N e , s per site and the parameters of the 
DFE. We obtained a low estimate for N e (posterior medians 350 for drug-treated and 
250 for control Influenza; Figure [t|A.) , which is expected given the bottleneck that the 
cells were subjected to in each transfer. While we obtained similar estimates for the y 
parameters for the drug-treated and for the control Influenza (posterior medians 0.44 and 
0.56, respectively), the a parameter was estimated to be much higher for the drug-treated 
than for the control Influenza (posterior medians 0.047 and 0.0071, respectively, Figure 
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|4^>). The resulting DFE was thus very different for the two conditions: the DFE for the 
drug-treated Influenza had a much heavier tail than the control (Figure [Ip). Posterior 
estimates for N e s per mutation also indicated that the drug-treated Influenza had more 
mutations under strong positive selection than the control (19% versus 3.5% of loci had 
P(N e s > 10) > 0.95, respectively; Figure [dp, Supplementary Figure S4). These results 
indicate that the drug treatment placed the Influenza population away from a fitness 
optimum, thus increasing the number of positively selected mutations with large effect 
sizes. Presumably these mutations confer resistance to the drug thus helping Influenza to 
reach a new fitness optimum. 

Our results for Influenza were qualitatively similar to those obtained by [26]. We 
obtained slightly larger estimates for N e (350 versus 226 for drug-treated and 250 versus 
176 for control Influenza). Our estimates for the parameters of the GPD were substantially 
different than [26] but resulted in qualitatively similar overall shapes of the DFE for both 
drug-treated and control experiments. These results underline the applicability of our 
method to a high-dimensional problem. In contrast to [26] who performed estimations in a 
3-step approach, combining a moment-based estimator for N e , ABC for s and a maximum 
likelihood approach for the GPD, our Bayesian framework allowed us to perform joint 
estimation and to obtain posterior distributions for all parameters in a single step. 


Conclusion 

Due to the difficulty to find analytically tractable likelihood solutions, statistical inference 
was often limited to models that made substantial approximations of reality. To address 
this problem, so-called likelihood-free approaches have been introduced that bypass the 
analytical evaluation of the likelihood function with computer simulations. While full- 
likelihood solutions generally have more power, likelihood-free methods have been used in 
many fields of science to overcome undesired model assumptions. 

Here we developed and implemented a novel likelihood-free, Markov chain Monte Carlo 
(MCMC) inference framework that scales naturally to high dimensions. This framework 
takes advantage of the observation that the information about one model parameter is 
often contained in a subset of the data, by integrating two key innovations: first, only a 
single parameter is updated at the time, and that update is accepted based on a subset 
of summary statistics sufficient for this parameter. We proved that this MCMC variant 
converges to the true joint posterior distribution under the standard assumptions. We 
further derived that for linear models, a one-dimensional function of summary statistics 
per parameter is sufficient. 

We then demonstrated the power of our framework through the application to multiple 
problems. First, our framework led to more accurate inference of the mean and standard 
deviation of a normal distribution than standard likelihood-free MCMC, suggesting that 
our framework is already competitive in models of low dimensionality. In high dimensions, 
the benefit was even more apparent. When applied to the problem of inferring parameters 
of a general linear model (GLM), for instance, we found our framework to be insensitive 
to the dimensionality, resulting in a performance similar to analytical solutions both in 
low and high dimensions. Finally, we used our framework to address the difficult and 
high-dimensional problem of inferring demography and selection jointly from genetic data. 
Specifically, and through simulations and an application to experimental data, we show 
that our framework enables the accurate joint estimation of the effective population size, 
the distribution of fitness effects of segregating mutations, as well as locus-specific selection 
coefficients from time series data. 
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Materials and methods 

Implementation 

We implemented the proposed ABC-PaSS framework into a new version of the software 
package ABCtoolbox m, which will be made available at the authors website and will be 
described elsewhere. 


Toy models: Normal distribution 

We performed simulations to assess the performance of ABC-MCMC and ABC-PaSS in esti¬ 
mating 9\ = /x and 62 = cr 2 for a univariate normal distribution. We used the sample mean 
x and sample variance S 2 of samples of size n as statistics. Recall that for non-informative 
priors the posterior distribution for /x is N(x, S 2 /n ) and the posterior distribution for cr 2 
is y 2 -distributed with n — 1 degrees of freedom. As /x and cr 2 are independent, we get the 
posterior density 


ir(/J, 1 (T~) fixjS^/n (/d 


n — 1 
S 2 


fx 2 ',ri— 1 


n — 1 

~g2 


-a 


In our simulations the sample size was n = 10 and the true parameters were given by 
/x = 0 and a 2 = 5. We performed 50 MCMC chains per simulation and chose effectively 
non-informative priors for fi t/[—10,10] and cr 2 ~ f7[0.1,15]. Our simulations were 
performed for a wide range of tolerances (from 0.1 to 0.9) and proposal ranges (from 0.1 
to 0.9). We did this exhaustive search in order to identify the combination of these tuning 
parameters that allows ABC-MCMC and ABC-PaSS to perform best in estimating /x and a 2 . We 
then recorded the minimum total variation distance (Li) between the true and estimated 
posteriors over these sets of tolerances and ranges and compared it between ABC-MCMC and 
ABC-PaSS. 


Toy models: GLM 


We considered linear models with m statistics s being a linear function of n 

0 : 


s = C 6 + e, e~A/’(0,/), 


m parameters 


where C is a square design matrix and the vector of errors e is multivariate normal. Under 
non-informative priors for the parameters 6 , their posterior distribution is multivariate 
normal 

6 \s ~ M (( C'C)- l C's , ( C'C )~ 1 ) . 

We set up the design matrices C in a cyclic manner to allow all statistics to have 
information on all parameters but their contributions to differ for each parameter, namely 
we set C = B ■ det(B ; B) 1 / 2n where 


( \/n 2/n 3/n 

n/n 1/n 2/n 

\ 2/10 3/10 4/10 


n/n \ 
n — 1/n 

2/n 

1/10 


The normalization factor in the definition of C was chosen such that the determinant 
of the posterior variance is constant and thus the widths of the marginal posteriors are 
comparable independently of the dimensionality n. We used all statistics for ABC-MCMC and 
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calculated a single linear combination of statistics per parameter for ABC-PaSS according 
to Theorem [2j For the estimation, we assumed that 6 = 0 and the priors are uniform 
t/[—100,100] for all parameters, which are effectively non-informative. We started the 
MCMC chains at a normal deviate N(6, 0.01/), i.e. around the true values of 6 . To ensure 
fair comparisons between methods, we performed simulations of 50 chains for a variety of 
tolerances (from 0.1 to 0.9) and proposal ranges (from 0.1 to 0.9) to choose the combination 
of these tuning parameters at which each method performed best. We run all our MCMC 
chains for 10 5 iterations per model parameter to account for model complexity. 

ABC-PaSS for estimating selection and demography 
Model 

Consider a vector £ of observed allele trajectories (sample allele frequencies) over l = 
1,..., L loci, as is commonly obtained in studies of experimental evolution. We assume 
these trajectories to be the result of both random drift as well as selection, parameterized 
by the effective population size N e and locus-specific selection coefficients s/, respectively, 
under the classic Wright-Fisher model with allelic fitnesses 1 and 1 + si. We further 
assume the locus-specific selection coefficients si follow a distribution of fitness effects 
(DFE) parameterized as a Generalized Pareto distribution (GPD) with mean y = 0, shape 
y and scale a. Our goal is thus to estimate the joint posterior distribution 

L 

n(N e ,si,...,s L ,x,<r\€) oc n [P(£z|iV e , si)n(si\x, u)]vr(lV e )7r(x)vr(cr) 

i=i 

To apply our ABC-PaSS framework to this problem, we approximate the likelihood term 
P(£;|lV e ,sz) numerically with simulations, while updating the hyper-parameters x and cr 
analytically. 


Summary statistics 

To summarize the data we used statistics originally proposed by Foil et al. [2TJ. Specif¬ 
ically, we first calculated for each locus individually a measure of the difference in allele 
frequency between consecutive time points as: 


, _ 1 Fs[ 1 - l/(2n)] - 2/n 
= t(l + Fs/A)[l-l/(n y )]' 


where 


Fs 


(x - y) 2 
*(!-*)’ 


x and y are the minor allele frequencies separated by t generations, 2 = (x + y )/2 and n 
is the harmonic mean of the sample sizes n x and n y . We then summed the Fs' values of 
all pairs of consecutive time points with increasing and decreasing allele frequencies into 
Fs'i and Fs'd, respectively [21]. Finally, we followed [28] and calculated boosted variants 
of the two statistics in order to take more complex relationships between parameters and 
statistics into account. The full set of statistics used per locus were Fi — {Fs'ii, Fs'di, 
Fs'if. Fs'd 2 , Fs'ii X Fs'di} ■ 

We next calculated parameter-specific linear combinations for N e and locus-specific si 
following the procedure developed above. To do so, we simulated allele trajectories of a 
single locus for different values of N e and s sampled from their prior. We then calculated 
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Fi for each simulation and performed a boxcox transformation to linearize the relationships 
between statistics and parameters [Sam. We then fit a linear model as outlined in Equation 
[2] to estimate the coefficients of an approximately sufficient linear combination of F for 
each parameter N e and s. This resulted in t s (Fi ) = / 3 S F\ and 7jv e (Fj) = Fi. To 
combine information across loci when updating N e . we then calculated 

L 

T N e (F) = ^ (3 Ne Fi, 

1=1 

where F = {Th, ... ,Fi}. In summary, we used the ABC approximation 

¥(^\N e ,Sj) «P(||r s (F,) ~t s (F 1o J\\ < 6 Sl ,\\r Ne (F)-T Ne (F obs )\\ < 5 N J\N e , Sj ). 

Simulations and Application 

We applied our framework to allele frequency data for the whole Influenza H1N1 genome 
obtained in a recently published evolutionary experiment |26| . In this experiment, In¬ 
fluenza A/Brisbane/59/2007 (H1N1) was serially amplified on Madin-Darby canine kidney 
(MDCK) cells for 12 passages of 72 hours each, corresponding to roughly 13 generations 
(doublings). After the three initial passages, samples were passed either in the absence of 
drug, or in the presence of increasing concentrations of the antiviral drug oseltamivir. At 
the end of each passage, samples were collected for whole genome high throughput pop¬ 
ulation sequencing. We obtained the raw data from http://bib.umassmed.edu/influenza/ 
and, following the original study [26], we downsampled it to 1000 haplotypes per time- 
point and filtered it to contain only loci for which sufficient data was available to calculate 
the Fs' statistics. Specifically, we included all loci with an allele frequency > 2% at > 2 
timepoints. There were 86 and 42 such loci for the control and drug-treated experiment, 
respectively. Further, we restricted our analysis of the data of the drug-treated experiment 
to the last nine time points during which drug was administered. 

We performed all our Wright-Fisher simulations with in-house C++ code implemented 
as a module of ABCtoolbox. We simulated 13 generations between timepoints and a sample 
of size 1000 per timepoint. We set the prior for N e uniform on the log 10 scale such 
that log 10 (A e ) ~ t/[1.5,4.5] and for the parameters of the GPD x ~ £7[—0.2,1] and for 
1 ogio(o‘) rsj U[— 2.5,—0.5]. For the simulations where no DFE was assumed, we set the 
prior of s ~ U[ 0,1]. 

As above, we run all our ABC-PaSS chains for 10 5 iterations per model parameter to 
account for model complexity. To ensure fast convergence, the ABC-PaSS implementa¬ 
tion benefited from an initial calibration step we originally developed for ABC-MCMC and 
implemented in ABCtoolbox [4]. Specifically, we first generated 10,000 simulations with 
values drawn randomly from the prior. For each parameter, we then selected the 1% sub¬ 
set of these simulations with the smallest distances to the observed data based on the 
linear-combination specific for that parameter. These accepted simulations were used to 
calibrate three important metrics prior to the MCMC run: first, we set the parameter- 
specific tolerances 5i to the largest distance among the accepted simulations. Second, we 
set the width of the parameter-specific proposal kernel to half of the standard deviation 
of the accepted parameter values. Third, we chose the starting value of the chain for each 
parameter as the accepted simulation with smallest distance. Each chain was then run 
for 1,000 iterations, and new starting values were chosen randomly among the accepted 
calibration simulations for those parameters for which no update was accepted. This was 
repeated until all parameters were updated at least once. 
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Appendix 

Proof for Theorem 1. The transition kernel IC(G,9') associated with the Markov chain 
is zero if 9 and 9' differ in more than one component. If 0_j = 0U for some index i, then 
we have 

K(0, O') = PiPi (9, G') + (1 - r(G))5 e (G') (3) 

where pi(9,9 r ) = qi(9'\9)F(T = ti :0bs \9')h(G,9') 1 5q is the Dirac mass in 9, and 

n 

r ( e ) = 

2—1 

We may assume without loss of generality that 

#W)< 1 

7r(9) qi (9'\9) ~ ' 

From 0 we conclude 

P(s = s obs \9 ) = P(Tj = t ifibs \9)gi{s obs ,9-i). 



Setting 

c := (J P(s = s obs \0)n{G)d9 
and keeping in mind that 9 i = 9'_ i and h{9\ 6) = 1, we get 



7r(0|s = s obs )pi{9,9') 


7r(0|s = s obs )q{6'\9)F(Ti = t ifibs \9')h(9, 9') 
c P(a = s o6s |0)7t(0M0'|0)P(T; = ti'Obs 

n(9)qi(9 \9) 

c P(T* = t ij0bs \9)gi(s 0bs , 9-i)F(Ti = t i)0bs \9')-K{9')qi[9\9') 
c P(Tj = ti^G^g^Sohs^'^FiTi = ti, obs \6)ir(9')qi(G\9 r ) 
c P(a = s obs \9')K{G')F(T L = t i>obs \6)K(9') qi {9\9')h(9\6) 
n(9'\s = s ohs )pi(9',9). 


From this and equation ([3]) follows readily that the transition kernel /C(-, •) satisfies the 
detailed balanced equation 


tt(0|s = s obs )lC(G, 9') = tt{G'\s = s obs )K{9',9) 


of the Metropolis-Hastings chain. 


□ 


Proof for Theorem 2. It is easy to check that the mean of t % is p, = t'^CG + c) and 
its variance is <r 2 = The covariance between s and r is given by 


S ST = E ((s — C9 — c)(n — pi)) 

= E ( ee'ti ) = X s ti. 

Consider the conditional multinormal distribution s|rj. Using the well-known formula for 
the variance and the mean of a conditional multivariate normal (see e.g. [29] . p. 63), we 
get that the covariance of s|rj is given by 

^s|r = — cr 2 X sr X(, T 
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and thus is independent of G. The mean of s\ri is 

H s \ T = CG + c + a^ 2 '£ ST t' i (s - CG - c). 
The part of this expression depending on 9i is 


(j _ 


Inserting ij = 1 Cj we obtain 


«f = (a ~ <*)«. = o- 

C i 2j s Cj / 

Thus the distribution of s|r* is independent of 9, and hence t* is sufficient for 6 t . 

To prove the second part of the theorem, we observe that r is given by the linear model 

r = C'Ej's = C'^CG + C'T,- l c + t] 

with rj = C'Ti~ l e. Using Cov(r j) = C ,/ X“ 1 C' we get for the posterior variance 


(C / S s - 1 (C'S7 1 C)- 1 C / E7 1 C + S,- 1 ) = (C'E^C + E, 

Similarly one sees that the posterior mean is Dd. 


-i\-i 


= D. 


□ 
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LI distance to true posterior (p) ^ 




Figure 1: Comparison of performance between ABC-MCMC (blue) and ABC-PaSS 
(red) in estimating the parameters of a normal distribution. (A,B): the average 
over 50 chains of the L\ distance between the true and estimated posterior distribution for 
/r (A) and a 2 (B) for different tolerances. The dashed horizontal line is the L\ distance 
between the prior and the true posterior distribution. (C,D): The estimated posterior 
distribution for /i (C) and cr 2 (D) using the tolerance that led to the minimum Li distance 
from the true posterior (black). The dashed vertical line indicates the true values of the 
parameters. 
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Figure 2: The performance of ABC-MCMC (blue) and ABC-PaSS (red) for a GLM 
model with different numbers of parameters. (A) the average L\ distance between 
the true and estimated posterior distribution for different tolerances. Solid and dashed lines 
are for a GLM with two and four parameters, respectively. (B) the minimum L\ distance 
from the true posterior over different tolerances for increasing numbers of parameters. 
(A,B) The dashed line is the L\ distance between the prior and the posterior distribution. 
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Figure 3: Accuracy in inferring demographic and selection parameters by 
ABC-PaSS. Shown are the true versus estimated posterior medians for parameters N e (A), 
s per locus (B), y and o of the Generalized Pareto distribution (C and D, respectively). 
Boxplots summarize results from 25 replicate simulations, each with 100 loci. Uniform 
priors over the whole ranges shown were used. (A, B): N e assumed in the simulations is 
represented as a color gradient of red (low N e ) to yellow (high N e ). (C,D): Parameters g 
and N e were fixed to 0 and 10 3 , respectively, logger was fixed to -1 (C) and y was fixed to 
0.5 (D). 
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Figure 4: Inferred demography and selection for experimental evolution of 

Infuenza. We show results for the no-drug (control) and drug-treated Influenza in grey 
and orange, respectively. Shown are the posterior distributions for logioN e (A) and logger 
and x (B). In panel C, we plotted the modal distribution of fitness effects (DFE) with thick 
lines by integrating over the posterior of its parameters. The thin lines represent the DFEs 
obtained by drawing 100 samples from the posterior of a and y. Dashed lines in panels A 
and C correspond to the prior distributions. In panel D, the posterior estimates for Nes per 
locus versus the position of the loci in the genome are shown. Open circles indicate non¬ 
significant loci whereas closed, thick circles indicate significant loci ( P(N e s > 10) > 0.95, 
dashed line). 


19 
























Supplementary Material 


20 



A « I 

• • • • 

% • •• • 

B 


-1-1- 

Z 1 

s i 

i* ••••» • 

w 

i i i 

so so- 


o - 

jR* r=0.90, P<2.2x10" 16 

- 

SjK*jfr? * r=0.80, P < 2.2 x 1 0 - 16 

_ 

n- 1 - 1 - 1 - 1 -r 

LO 

h-1-1-1-1-1-H 


0.0 0.2 0.4 0.6 0.8 1.0 ' 1.5 2.5 3.5 4.5 

s logi 0 N e 


Figure SI: Relationship between parameters and linear combinations. The re¬ 
lationship between the parameters s and log\oN e with respective linear combinations of 
statistics t s and tn b calculated for a set of 10 4 simulations and assuming priors used for 
the Influenza application presented in the main text. The Pearson correlation coefficient 
(?’) and the corresponding P-value are shown in each panel. 
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Figure S2: Directed acyclic graph describing the Wright-Fisher model examined 
in this study. Solid circles represent parameters to be estimated. The dashed square 
represents the full data, which is summerized here by a vector of statistics F[, indicated 
by a solid square. Nodes contained in the plate are repeated for each locus l £ {1,..., L} 
times. 
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Figure S3: Accuracy in estimating N e and DFE parameters a and y jointly. 

(A,B,C) Sets of simulations of 100 loci were conducted for combinations of parameters a 
and x over a grid from their prior range and we evaluated the median approximation error 
(e=estimate-true) over 25 replicates. Color gradients indicate the extent of overestimation 
(red) or underestimation (blue) of each parameter. These results suggest very high accuracy 
when estimating N e with maximum e ~ 0.03 or 1% of the prior range and rather low for 
a (about 10% of the prior range). In contrast, e is rather large for y, spanning up to 75% 
of the prior range. This is due to several combinations of y and a leading to very similar 
shapes of the truncated GPD. This is illustrated in panels D, E anf F, where we show 
the true (dashed black line) versus estimated (red) DFE obtained for 25 replicates using 
parameter combinations of chi and a as indicated in panels B and C. 
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Figure S4: Allele trajectories and posterior estimates for N e s for control and 
drug-treated Influenza. Non-significant loci are colored grey and significant loci are 
colored with a unique color for each locus. 


24 


































































